Imports System.Drawing
Imports System.Text
Imports DotSpatial.Topology
Imports MWPoint = DotSpatial.Topology.Point

Friend Class TileCalculator
    ' Adapted from methods at http://msdn.microsoft.com/en-us/library/bb259689.aspx 


    Public Const MinLatitude As Double = -85.05112878
    Public Const MaxLatitude As Double = 85.05112878
    Public Const MinLongitude As Double = -180
    Public Const MaxLongitude As Double = 180
    Public Const EarthRadiusKms As Double = 6378137

    Public Const MaxWebMercX As Double = 20037497.2108402
    Public Const MinWebMercX As Double = -20037497.2108402
    Public Const MaxWebMercY As Double = 20037508.3430388
    Public Const MinWebMercY As Double = -20037508.3430388

    ''' <summary>
    '''
    ''' </summary>
    ''' <param name="envelope"></param>
    ''' <param name="rectangle"></param>
    ''' <returns></returns>
    Public Shared Function DetermineZoomLevel(envelope As Envelope, rectangle As Rectangle) As Integer
        Dim metersAcross As Double = EarthRadiusKms * envelope.Width * Math.PI / 180
        'find the arc length represented by the displayed map
        metersAcross *= Math.Cos(envelope.Center().Y * Math.PI / 180)
        'correct for the center latitude
        Dim metersAcrossPerPixel As Double = metersAcross / rectangle.Width
        'find the resolution in meters per pixel
        'find zoomlevel such that metersAcrossPerPix is close
        For i As Integer = 2 To 18
            Dim groundRes As Double = GroundResolution(envelope.Center().Y, i)

            If metersAcrossPerPixel > groundRes Then
                'System.Diagnostics.Debug.WriteLine("metersPerPixel: " + metersAcrossPerPixel);
                'System.Diagnostics.Debug.WriteLine("groundRes:      " + groundRes);
                'double ratio = metersAcrossPerPixel / groundRes;
                'System.Diagnostics.Debug.WriteLine("ratio: " + ratio);

                'fix zoom level..
                'changed to a slightly lower zoom level to increase readability
                If i > 2 AndAlso i < 18 Then
                    Return i - 1
                End If

                'MessageBox.Show("MAPP: "+metersAcrossPerPixel+" , zoom: "+(i-1));
                Return i
            End If
        Next

        Return -1
    End Function

    ''' <summary>
    ''' Clips a number to the specified minimum and maximum values.
    ''' </summary>
    ''' <param name="n">The number to clip.</param>
    ''' <param name="minValue">Minimum allowable value.</param>
    ''' <param name="maxValue">Maximum allowable value.</param>
    ''' <returns>The clipped value.</returns>
    Public Shared Function Clip(n As Double, minValue As Double, maxValue As Double) As Double
        Return Math.Min(Math.Max(n, minValue), maxValue)
    End Function

    ''' <summary>
    ''' Determines the map width and height (in pixels) at a specified level
    ''' of detail.
    ''' </summary>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    ''' <returns>The map width and height in pixels.</returns>
    Public Shared Function MapSize(levelOfDetail As Integer) As UInteger
        Return CUInt(256) << levelOfDetail
    End Function

    ''' <summary>
    ''' Determines the ground resolution (in meters per pixel) at a specified
    ''' latitude and level of detail.
    ''' </summary>
    ''' <param name="latitude">Latitude (in degrees) at which to measure the
    ''' ground resolution.  </param>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    ''' <returns>The ground resolution, in meters per pixel.</returns>
    Public Shared Function GroundResolution(latitude As Double, levelOfDetail As Integer) As Double
        latitude = Clip(latitude, MinLatitude, MaxLatitude)
        Return Math.Cos(latitude * Math.PI / 180) * 2 * Math.PI * EarthRadiusKms / MapSize(levelOfDetail)
    End Function

    ''' <summary>
    ''' Determines the map scale at a specified latitude, level of detail,
    ''' and screen resolution.
    ''' </summary>
    ''' <param name="latitude">Latitude (in degrees) at which to measure the
    ''' map scale.  </param>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    ''' <param name="screenDpi">Resolution of the screen, in dots per inch.</param>
    ''' <returns>The map scale, expressed as the denominator N of the ratio 1 : N.</returns>
    Public Shared Function MapScale(latitude As Double, levelOfDetail As Integer, screenDpi As Integer) As Double
        Return GroundResolution(latitude, levelOfDetail) * screenDpi / 0.0254
    End Function

    ''' <summary>
    ''' Converts a point from latitude/longitude WGS-84 coordinates (in degrees)
    ''' into pixel XY coordinates at a specified level of detail.
    ''' </summary>
    ''' <param name="latitude">Latitude of the point, in degrees.</param>
    ''' <param name="longitude">Longitude of the point, in degrees.</param>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    Public Shared Function LatLongToPixelXY(latitude As Double, longitude As Double, levelOfDetail As Integer) As MWPoint
        latitude = Clip(latitude, MinLatitude, MaxLatitude)
        longitude = Clip(longitude, MinLongitude, MaxLongitude)

        Dim x = (longitude + 180) / 360
        Dim sinLatitude = Math.Sin(latitude * Math.PI / 180)
        Dim y = 0.5 - Math.Log((1 + sinLatitude) / (1 - sinLatitude)) / (4 * Math.PI)

        Dim mapSize__1 = MapSize(levelOfDetail)

        Dim pixelXY = New MWPoint(CInt(Math.Truncate(Clip(x * mapSize__1 + 0.5, 0, mapSize__1 - 1))), CInt(Math.Truncate(Clip(y * mapSize__1 + 0.5, 0, mapSize__1 - 1))))

        Return pixelXY
    End Function

    ''' <summary>
    ''' Converts a point from latitude/longitude WGS-84 coordinates (in degrees)
    ''' into pixel XY coordinates at a specified level of detail.
    ''' </summary>
    ''' <param name="coord">Latitude, Longitude Coordinate of the point, in degrees.</param>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    Public Shared Function LatLongToPixelXY(coord As Coordinate, levelOfDetail As Integer) As MWPoint
        Return LatLongToPixelXY(coord.Y, coord.X, levelOfDetail)
    End Function

    ''' <summary>
    ''' Converts a pixel from pixel XY coordinates at a specified level of detail
    ''' into latitude/longitude WGS-84 coordinates (in degrees).
    ''' </summary>
    ''' <param name="pixelX">X coordinate of the point, in pixels.</param>
    ''' <param name="pixelY">Y coordinates of the point, in pixels.</param>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    Public Shared Function PixelXYToLatLong(pixelX As Integer, pixelY As Integer, levelOfDetail As Integer) As Coordinate
        Dim mapSize__1 As Double = MapSize(levelOfDetail)
        Dim x As Double = (Clip(pixelX, 0, mapSize__1 - 1) / mapSize__1) - 0.5
        Dim y As Double = 0.5 - (Clip(pixelY, 0, mapSize__1 - 1) / mapSize__1)

        Dim latitude As Double = 90 - 360 * Math.Atan(Math.Exp(-y * 2 * Math.PI)) / Math.PI
        Dim longitude As Double = 360 * x

        Return New Coordinate(longitude, latitude)
    End Function

    ''' <summary>
    ''' Converts pixel XY coordinates into tile XY coordinates of the tile containing
    ''' the specified pixel.
    ''' </summary>
    ''' <param name="pixelX">Pixel X coordinate.</param>
    ''' <param name="pixelY">Pixel Y coordinate.</param>
    Public Shared Function PixelXYToTileXY(pixelX As Integer, pixelY As Integer) As MWPoint
        Return New MWPoint(pixelX \ 256, pixelY \ 256)
    End Function

    ''' <summary>
    ''' Converts pixel XY coordinates into tile XY coordinates of the tile containing
    ''' the specified pixel.
    ''' </summary>
    ''' <param name="point">Pixel X,Y point.</param>
    Public Shared Function PixelXYToTileXY(point As MWPoint) As MWPoint
        Return PixelXYToTileXY(CInt(point.X), CInt(point.Y))
    End Function

    ''' <summary>
    ''' Converts tile XY coordinates into pixel XY coordinates of the upper-left pixel
    ''' of the specified tile.
    ''' </summary>
    ''' <param name="tileX">Tile X coordinate.</param>
    ''' <param name="tileY">Tile Y coordinate.</param>
    Public Shared Function TileXYToTopLeftPixelXY(tileX As Integer, tileY As Integer) As MWPoint
        Return New MWPoint(tileX * 256, tileY * 256)
    End Function

    ''' <summary>
    ''' Converts tile XY coordinates into pixel XY coordinates of the bottom-right pixel
    ''' of the specified tile.
    ''' </summary>
    ''' <param name="tileX"></param>
    ''' <param name="tileY"></param>
    ''' <returns></returns>
    Public Shared Function TileXYToBottomRightPixelXY(tileX As Integer, tileY As Integer) As MWPoint
        Return New MWPoint((tileX + 1) * 256 - 1, (tileY + 1) * 256 - 1)
    End Function

    ''' <summary>
    ''' Converts a WGS-84 Lat/Long coordinate to the tile XY of the tile containing
    ''' that point at the given levelOfDetail
    ''' </summary>
    ''' <param name="coord">WGS-84 Lat/Long</param>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    ''' <returns>Tile XY Point</returns>
    Public Shared Function LatLongToTileXY(coord As Coordinate, levelOfDetail As Integer) As MWPoint
        Dim pixelXY As MWPoint = LatLongToPixelXY(coord.Y, coord.X, levelOfDetail)
        Dim tileXY As MWPoint = PixelXYToTileXY(pixelXY)

        Return tileXY
    End Function

    ''' <summary>
    ''' Converts tile XY coordinates into a QuadKey at a specified level of detail.
    ''' </summary>
    ''' <param name="tileX">Tile X coordinate.</param>
    ''' <param name="tileY">Tile Y coordinate.</param>
    ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
    ''' to 23 (highest detail).  </param>
    ''' <returns>A string containing the Bing QuadKey.</returns>
    Public Shared Function TileXYToBingQuadKey(tileX As Integer, tileY As Integer, levelOfDetail As Integer) As String
        Dim quadKey = New StringBuilder()
        For i As Integer = levelOfDetail To 1 Step -1
            Dim digit As Char = "0"
            Dim mask = 1 << (i - 1)
            If (tileX And mask) <> 0 Then
                digit = Chr(Asc(digit) + 1)
            End If
            If (tileY And mask) <> 0 Then
                digit = Chr(Asc(digit) + 2)
            End If
            quadKey.Append(digit)
        Next
        Return quadKey.ToString()
    End Function

    ''' <summary>
    ''' Converts a Bing QuadKey into tile XY coordinates.
    ''' </summary>
    ''' <param name="quadKey">QuadKey of the tile.</param>
    ''' <param name="tileX">Output parameter receiving the tile X coordinate.</param>
    ''' <param name="tileY">Output parameter receiving the tile Y coordinate.</param>
    ''' <param name="levelOfDetail">Output parameter receiving the level of detail.</param>
    Public Shared Sub BingQuadKeyToTileXY(quadKey As String, ByRef tileX As Integer, ByRef tileY As Integer, ByRef levelOfDetail As Integer)
        tileX = InlineAssignHelper(tileY, 0)
        levelOfDetail = quadKey.Length
        For i As Integer = levelOfDetail To 1 Step -1
            Dim mask As Integer = 1 << (i - 1)
            Select Case quadKey(levelOfDetail - i)
                Case "0"c
                    Exit Select

                Case "1"c
                    tileX = tileX Or mask
                    Exit Select

                Case "2"c
                    tileY = tileY Or mask
                    Exit Select

                Case "3"c
                    tileX = tileX Or mask
                    tileY = tileY Or mask
                    Exit Select
                Case Else

                    Throw New ArgumentException("Invalid QuadKey digit sequence.")
            End Select
        Next
    End Sub

    ''' <summary>
    ''' Takes a 2-dimensional array of tiles and stitches it into a single Bitmap for display on the Map
    ''' </summary>
    ''' <param name="tiles">2-dimensional array of tiles, [x by y]</param>
    ''' <returns>Bitmap of the tiles stitched together</returns>
    Public Shared Function StitchTiles(tiles As Tile(,)) As Bitmap
        Dim width = tiles.GetLength(0) * 256
        Dim height = tiles.GetLength(1) * 256

        'create a bitmap to hold the combined image
        Dim finalImage = New Bitmap(width, height)

        'get a graphics object from the image so we can draw on it
        Using g As Graphics = Graphics.FromImage(finalImage)
            'set background color
            g.Clear(Color.Black)

            'go through each image and "draw" it on the final image

            For y As Integer = 0 To tiles.GetLength(1) - 1
                For x As Integer = 0 To tiles.GetLength(0) - 1
                    If tiles(x, y) IsNot Nothing Then
                        Dim tile = tiles(x, y).Bitmap

                        g.DrawImage(tile, New Rectangle(x * 256, y * 256, tile.Width, tile.Height))
                    End If
                Next
            Next
        End Using

        Return finalImage
    End Function
    Private Shared Function InlineAssignHelper(Of T)(ByRef target As T, value As T) As T
        target = value
        Return value
    End Function
End Class
